/*********************************************************************************************************
 *  ------------------------------------------------------------------------------------------------------
 *  file description
 *  ------------------------------------------------------------------------------------------------------
 *         \file  vmath.c
 *         \unit  vmath
 *        \brief  Similar to the C standard library math
 *       \author  Lamdonn
 *      \version  v1.0.0
 *      \license  GPL-2.0
 *    \copyright  Copyright (C) 2023 Lamdonn.
 ********************************************************************************************************/
#include "vmath.h"

static double _ln[100] = { 0.0,0.009950330853168,0.019802627296180,0.029558802241544,0.039220713153281,0.048790164169432,0.058268908123976,0.067658648473815,0.076961041136128,0.086177696241052,0.095310179804325,0.104360015324243,0.113328685307003,0.122217632724249,0.131028262406404,0.139761942375159,0.148420005118273,0.157003748809665,0.165514438477574,0.173953307123438,0.182321556793955,0.190620359608650,0.198850858745165,0.207014169384326,0.215111379616946,0.223143551314210,0.231111720963387,0.239016900470500,0.246860077931526,0.254642218373581,0.262364264467491,0.270027137213060,0.277631736598280,0.285178942233663,0.292669613962820,0.300104592450338,0.307484699747961,0.314810739840034,0.322083499169114,0.329303747142601,0.336472236621213,0.343589704390077,0.350656871613170,0.357674444271816,0.364643113587910,0.371563556432483,0.378436435720245,0.385262400790645,0.392042087776024,0.398776119957368,0.405465108108165,0.412109650826833,0.418710334858185,0.425267735404344,0.431782416425538,0.438254930931156,0.444685821261446,0.451075619360217,0.457424847038876,0.463734016232140,0.470003629245736,0.476234178996372,0.482426149244293,0.488580014818671,0.494696241836107,0.500775287912490,0.506817602368452,0.512823626428664,0.518793793415168,0.524728528934982,0.530628251062171,0.536493370514569,0.542324290825362,0.548121408509688,0.553885113226438,0.559615787935423,0.565313809050061,0.570979546585738,0.576613364303994,0.582215619852664,0.587786664902119,0.593326845277735,0.598836501088704,0.604315966853330,0.609765571620895,0.615185639090234,0.620576487725110,0.625938430866496,0.631271776841858,0.636576829071551,0.641853886172395,0.647103242058539,0.652325186039691,0.657520002916795,0.662687973075237,0.667829372575656,0.672944473242426,0.678033542749898,0.683096844706444,0.688134638736401 };

static double _exp1[71] = { 1.0,2.202646579480672e+004,4.851651954097903e+008,1.068647458152446e+013,2.353852668370200e+017,5.184705528587072e+021,1.142007389815684e+026,2.515438670919167e+030,5.540622384393510e+034,1.220403294317841e+039,2.688117141816136e+043,5.920972027664670e+047,1.304180878393632e+052,2.872649550817832e+056,6.327431707155585e+060,1.393709580666380e+065,3.069849640644242e+069,6.761793810485009e+073,1.489384200781838e+078,3.280587015384671e+082,7.225973768125749e+086,1.591626640377924e+091,3.505790975238748e+095,7.722018499983836e+099,1.700887763567586e+104,3.746454614502673e+108,8.252115441813891e+112,1.817649385139100e+117,4.003639200871785e+121,8.818602191274965e+125,1.942426395241256e+130,4.278478855371123e+134,9.423976816163585e+138,2.075769029922787e+143,4.572185553551339e+147,1.007090887028080e+152,2.218265297538556e+156,4.886054470003974e+160,1.076225116551050e+165,2.370543571722357e+169,5.221469689764144e+173,1.150105235202100e+178,2.533275362360718e+182,5.579910311786494e+186,1.229057036206545e+191,2.707178276786998e+195,5.962956971409261e+199,1.313428677666503e+204,2.893019184253945e+208,6.372298810568915e+212,1.403592217852838e+217,3.091617597639242e+221,6.809740926502327e+225,1.499945255890989e+230,3.303849287296548e+234,7.277212331783397e+238,1.602912685075726e+243,3.530650142988227e+247,7.776774460795963e+251,1.712948566546487e+256,3.773020300929940e+260,8.310630260154467e+264,1.830538131585780e+269,4.032028554146358e+273,8.881133903158874e+277,1.956199921370272e+282,4.308817065586588e+286,9.490801171122244e+290,2.090488073610356e+295,4.604606404782990e+299 };

static double _exp2[70] = { 1.0,4.539992976248485e-005,2.061153622438558e-009,9.357622968840175e-014,4.248354255291589e-018,1.928749847963918e-022,8.756510762696520e-027,3.975449735908647e-031,1.804851387845415e-035,8.194012623990515e-040,3.720075976020836e-044,1.688911880224532e-048,7.667648073722000e-053,3.481106839904311e-057,1.580420060273613e-061,7.175095973164411e-066,3.257488532207521e-070,1.478897505643213e-074,6.714184288211594e-079,3.048234950971857e-083,1.383896526736738e-087,6.282880511239462e-092,2.852423339163565e-096,1.294998192508984e-100,5.879282698245269e-105,2.669190215541276e-109,1.211810483082857e-113,5.501611081740457e-118,2.497727566915251e-122,1.133966561037746e-126,5.148200222412014e-131,2.337279285007143e-135,1.061123153746351e-139,4.817491664943076e-144,2.187137832197718e-148,9.929590396264980e-153,4.508027065606742e-157,2.046641121459268e-161,9.291736316326398e-166,4.218441761327482e-170,1.915169596714006e-174,8.694856517406230e-179,3.947458751851265e-183,1.792143500743535e-187,8.136318905805023e-192,3.693883068487256e-196,1.677020318601535e-200,7.613660467476964e-205,3.456596504588617e-209,1.569292385255739e-213,7.124576406741286e-218,3.234552684535111e-222,1.468484646909508e-226,6.666909982697905e-231,3.026772449472940e-235,1.374152566130957e-239,6.238642998528377e-244,2.832339539464062e-248,1.285880161551771e-252,5.837886901742308e-257,2.650396553004311e-261,1.203278173491277e-265,5.462874456123503e-270,2.480141166092796e-274,1.125982347416602e-278,5.111951948651156e-283,2.320822594179601e-287,1.053651827669418e-291,4.783571897030535e-296,2.171738281389827e-300 };

static double _exp3[101] = { 1.0,1.105170918075648,1.221402758160170,1.349858807576003,1.491824697641270,1.648721270700128,1.822118800390509,2.013752707470477,2.225540928492467,2.459603111156949,2.718281828459045,3.004166023946433,3.320116922736547,3.669296667619244,4.055199966844675,4.481689070338065,4.953032424395117,5.473947391727202,6.049647464412949,6.685894442279273,7.389056098930653,8.166169912567655,9.025013499434127,9.974182454814727,11.023176380641610,12.182493960703484,13.463738035001704,14.879731724872849,16.444646771097069,18.174145369443085,20.085536923187693,22.197951281441664,24.532530197109384,27.112638920657929,29.964100047397064,33.115451958692375,36.598234443678059,40.447304360067470,44.701184493300914,49.402449105530280,54.598150033144336,60.340287597362057,66.686331040925211,73.699793699595844,81.450868664968141,90.017131300521811,99.484315641933776,109.94717245212343,121.51041751873476,134.28977968493530,148.41315910257634,164.02190729990139,181.27224187515074,200.33680997479112,221.40641620418637,244.69193226421953,270.42640742615157,298.86740096705898,330.29955990964714,365.03746786532696,403.42879349273295,445.85777008251438,492.74904109325325,544.57191012592557,601.84503787207802,665.14163304435715,735.09518924196743,812.40582516753682,897.84729165041040,992.27471560501738,1096.6331584284490,1211.9670744925654,1339.4307643944051,1480.2999275845305,1635.9844299959097,1808.0424144560438,1998.1958951040961,2208.3479918871835,2440.6019776244702,2697.2823282684762,2980.9579870416910,3294.4680752837994,3640.9503073323067,4023.8723938222556,4447.0667476997942,4914.7688402990643,5431.6595913629008,6002.9122172609323,6634.2440062777841,7331.9735391558779,8103.0839275752542,8955.2927034823661,9897.1290587437506,10938.019208164997,12088.380730216773,13359.726829661635,14764.781565577005,16317.607198015130,18033.744927828171,19930.370438229907,22026.465794806718 };

static double _exp4[101] = { 1.0,1.001000500166708,1.002002001334000,1.003004504503377,1.004008010677342,1.005012520859401,1.006018036054065,1.007024557266849,1.008032085504274,1.009040621773868,1.010050167084168,1.011060722444720,1.012072288866078,1.013084867359809,1.014098458938492,1.015113064615719,1.016128685406095,1.017145322325241,1.018162976389794,1.019181648617408,1.020201340026756,1.021222051637529,1.022243784470438,1.023266539547218,1.024290317890622,1.025315120524429,1.026340948473442,1.027367802763489,1.028395684421425,1.029424594475131,1.030454533953517,1.031485503886523,1.032517505305119,1.033550539241306,1.034584606728118,1.035619708799623,1.036655846490924,1.037693020838157,1.038731232878498,1.039770483650158,1.040810774192388,1.041852105545480,1.042894478750763,1.043937894850613,1.044982354888444,1.046027859908717,1.047074410956937,1.048122009079656,1.049170655324471,1.050220350740028,1.051271096376024,1.052322893283204,1.053375742513365,1.054429645119356,1.055484602155080,1.056540614675494,1.057597683736611,1.058655810395500,1.059714995710288,1.060775240740159,1.061836546545360,1.062898914187195,1.063962344728034,1.065026839231306,1.066092398761505,1.067159024384193,1.068226717165993,1.069295478174600,1.070365308478774,1.071436209148346,1.072508181254217,1.073581225868358,1.074655344063814,1.075730536914703,1.076806805496220,1.077884150884632,1.078962574157284,1.080042076392600,1.081122658670083,1.082204322070315,1.083287067674959,1.084370896566760,1.085455809829549,1.086541808548238,1.087628893808826,1.088717066698399,1.089806328305129,1.090896679718278,1.091988122028198,1.093080656326330,1.094174283705210,1.095269005258466,1.096364822080817,1.097461735268082,1.098559745917174,1.099658855126103,1.100759063993979,1.101860373621011,1.102962785108508,1.104066299558882,1.105170918075648 };

double v_fmod(double x, double y)
{
    return x - ((int)(x / y)) * y;
}

static int factorial(int n)
{
    int res = 1, i;
    for (i = 1; i <= n; i++)
    {
        res *= i;
    }
    return res;
}

static double sector(double x)
{
    int i;
    double sum = 0, t;
    for (i = 1; i <= 10; i += 2)
    {
        t = v_pow(x, i) / factorial(i);
        if (i % 4 == 1) sum += t;
        else if (i % 4 == 3) sum -= t;
    }
    return sum;
}

double v_sin(double x)
{
    x = v_fmod(x, 2 * PI);
    if (x < 0) return -v_sin(-x); // sin(x) = -sin(-x)
    else if (x == 0) return 0;
    else if (x < PI * 0.5) return sector(x);
    else if (x == PI * 0.5) return 1;
    else if (x <= PI) return v_sin(PI - x); // sin(x) = sin(pi - x)
    else if (x <= PI * 2) return -v_sin(2 * PI - x); // f(x) = -sin(2pi - x)
    return 0;
}

double v_cos(double x)
{
    x = v_fmod(x, 2 * PI);
    return v_sin(x + PI * 0.5); // cos(x) = sin(x + 0.5pi)
}

double v_sec(double x)
{
    return 1.0 / v_cos(x);
}

double v_csc(double x)
{
    return 1.0 / v_sin(x);
}

double v_tan(double x)
{
    double s = v_sin(x);
    double c = v_cos(x);
    if (c == 0) return s > 0 ? INFINITY : -INFINITY;
    return s / c;
}

double v_cot(double x)
{
    return v_tan(0.5 * PI - x);
}

static inline double i_ln(double x)
{
    long long* b = (long long*)(&x);
    int c = ((*b) >> 52) - 1023;
    double rst = ((double)c) * 0.69314718055994530941723212145818;
    int d;
    double e;

    *b = ((*b) & 0x000fffffffffffff) | 0x3ff0000000000000;
    d = (int)(x * 100.0);
    rst += _ln[d - 100];
    x = (x * 100.0 / (double)d) - 1.0;
    e = x * x;

    return (rst + x - e / 2.0 + e * x / 3.0 - e * e / 4.0);
}

double v_ln(double x)
{
    if (x <= 0.0) return -INFINITY;
    return i_ln(x);
}

#if 0
#define N 100
double m_ln(double x) {
    double res = 0.0;
    double h = (x - 1) / N; // N为分段数，可以根据需要调整
    for (int i = 1; i <= N; i++) {
        double ti = 1 + i * h;
        res += h * (1 / ti);
    }
    return res;
}
#endif

double v_lg(double x)
{
    if (x <= 0.0) return -INFINITY;
    return i_ln(x) / 2.3025850929940456840179914546844;
}

double v_log2(double x)
{
    if (x <= 0.0) return -INFINITY;
    return i_ln(x) / 0.69314718055994530941723212145818;
}

double v_loga(double x, double y)
{
    if ((x == 0.0) || (y <= 0.0)) return -INFINITY;
    return i_ln(y) / i_ln(x);
}

double v_exp(double x)
{
    if ((x > 700.0) || (x < -700)) return -INFINITY;
    double f = x;
    int b = (int)(f * 1000.0);
    f -= (double)(b) / 1000.0;
    int c = b / 10000 - (x < 0); b -= c * 10000;
    double rst = (x < 0) ? _exp2[0 - c] : _exp1[c];

    int d = b / 100; int e = b % 100;
    rst = rst * _exp3[d] * _exp4[e];

    double g = f * f;
    return (1 + f + g / 2.0 + g * f / 6.0 + g * g / 24.0) * rst;
}

double v_pow(double x, double y)
{
    if (x <= 0.0) return -INFINITY;
    return v_exp(y * i_ln(x));
}

double v_sigmoid(double x)
{
    return (1.0 / (1.0 + v_exp(0.0 - x)));
}

double v_dsigmoid(double x)
{
    double s = (1.0 / (1.0 + v_exp(0.0 - x)));
    return s * (1 - s);
}

double v_sinh(double x)
{
    double a = v_exp(x);
    return (a - 1.0 / a) / 2.0;
}

double v_cosh(double x)
{
    double a = v_exp(x);
    return (a + 1.0 / a) / 2.0;
}

double v_tanh(double x)
{
    double a = v_exp(x); double b = 1.0 / a;
    return (a - b) / (a + b);
}

double v_coth(double x)
{
    double a = v_exp(x); double b = 1.0 / a;
    return (a + b) / (a - b);
}

double v_sech(double x)
{
    double a = v_exp(x);
    return 2.0 / (a + 1.0 / a);
}

double v_csch(double x)
{
    double a = v_exp(x);
    return 2.0 / (a - 1.0 / a);
}

double v_sqrt(double x)
{
    if (x < 0.0) return -INFINITY;
    return v_pow(x, 0.5);
}

double v_gaussian(double x, double mean, double variance)
{
    double b = (x - mean); b = 0 - b * b / (variance + variance);
    return v_exp(b) / (v_sqrt(variance) * 2.506628274631000502415765284811);
}
